This script analyzes bird community data collected for the Woodland Restoration Survey in 2022. Goals of the study are to assess if (1) bird communities in restored woodlots are progressing towards mature forest reference sites., (2) if actively restored sites have different trajectories than passively restored sites, and (3) if there are any differences in how birds respond to different woodland restoration methods (large bareroot plantings vs. seedling plantings).
Helpful references used to guide this analysis:
library(tidyverse)
library(sf)
library(vegan)
library(goeveg) # needed for dimcheckMDS to create scree plot.
library(here)
library(pairwiseAdonis)
library(RVAideMemoire)
`%!in%` <- Negate(`%in%`)
source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
# source(here("Scripts", "4_Covariate_Data_Prep.R"))
# source(here("Scripts", "5_Prep_Species.R"))wrs <- wrs %>%
filter(Temp_issue != 1 | is.na(Temp_issue)) %>% # remove rows that need to be fixed
filter(Species != "") %>% # remove blank species (temporary)
filter(Flyover != 1 | is.na(Flyover)) %>% # remove flyovers
filter(!Initial_detection_time_interval %in% c("After count", "Before count")) %>% # remove detection occurring before/after count
filter(!Initial_detection_distance_bin %in% c("Beyond woodlot"))
# Impute missing distance bins from best guess
wrs <- wrs %>%
mutate(Initial_detection_distance_bin = if_else(Initial_detection_distance_bin == "", ">50 m but in woodlot", Initial_detection_distance_bin))
# Split species code and common name
wrs <- wrs %>%
separate(Species, c("Sp", "Common_name"), sep = ": ")
# Drop unidentified birds
wrs <- wrs[!grepl("Unidentified", wrs$Common_name), ]
# Get unique sites (create standalone df)
wrs_sites <- unique(wrs$Site_ID_full) %>%
as.data.frame() %>%
rename("Site_ID_full" = ".")
# Get unique surveys (add to main df and create standalone)
wrs <- wrs %>%
mutate(Survey_ID_full = paste0(Site_ID_full, " - ", Date))
wrs_surveys <- unique(wrs$Survey_ID_full) %>%
as.data.frame() %>%
rename("Survey_ID_full" = ".") %>%
separate(col = Survey_ID_full, remove = FALSE, sep = " - ", into = c("SurveyType", "Park", "Site", "Date")) %>%
mutate(Site_ID_full = paste(SurveyType, Park, Site, sep = " - ")) %>%
dplyr::select(Survey_ID_full, Site_ID_full)
# TODO add survey covariates to survey table (wind, time, etc)
# Join site info from shapefile
# Make site IDs in site info df match those in bird data
pts <- st_read(here("Projects", "Woodland_Bird_Surveys", "Survey_locations_attributed.shp")) %>%
mutate(Park = ifelse(Park == "Cleary Lake", "Cleary", Park))
site_info <- pts
# Drop the truncated park names from part after final "-"
site_info$Site_ID_full <- gsub("CH", "", site_info$SurveyID)
site_info$Site_ID_full <- gsub("BK", "", site_info$Site_ID_full)
site_info$Site_ID_full <- gsub("LR", "", site_info$Site_ID_full)
# Remove non-trailing zeros
site_info$Site_ID_full <- gsub("(?<![0-9])0+", "", site_info$Site_ID_full, perl = TRUE)
# Tweak survey type names in full ID to match
site_info$Site_ID_full <- gsub("WRS Point Count", "Woodland Restoration Survey", site_info$Site_ID_full)
site_info$Site_ID_full <- gsub("Point Count MB", "Mountain Bike Trail Survey", site_info$Site_ID_full)
site_info <- site_info %>%
dplyr::select(Site_ID_full, Type, YearPlant)
# Manually append info on points added later
# TODO add these to actual shapefile and confirm values
new_site_data <- data.frame(
Site_ID_full = c("Woodland Point Count - Baker - 1", "Woodland Point Count - Baker - 2"),
Type = c("Mature forest", "Mature forest"),
YearPlant = c(1953, 1953)
)
site_info <- bind_rows(site_info, new_site_data)
wrs_sites <- left_join(wrs_sites, site_info, by = join_by(Site_ID_full))
# Arrange data
wrs_surveys <- left_join(wrs_surveys, site_info, by = join_by(Site_ID_full)) %>%
arrange(Survey_ID_full) %>%
dplyr::select(-c(Survey_ID_full, Site_ID_full, geometry)) %>%
mutate(
Time = 2022 - YearPlant,
Type = as.factor(Type)
) %>%
dplyr::select(Time, Type)
wrs_surveys$Type <- factor(wrs_surveys$Type, levels = c("Mature forest", "LBR", "Seedling", "Natural succession")) # reorder factor levels, moving mature first
wrs_species <- wrs %>%
group_by(Survey_ID_full, Sp) %>%
summarize(Count = n()) %>%
spread(Sp, Count, fill = 0) %>%
arrange(Survey_ID_full) %>%
ungroup() %>%
dplyr::select(-Survey_ID_full)sites_for_map <- wrs_sites %>%
separate(Site_ID_full, sep = " - ", into = c("Survey_type", "Park", "ExistingID")) %>%
st_as_sf()
source(here("Scripts", "986_Mapping_Functions.R"))
forestry_sites <- st_read(here("Data", "Geospatial", "Vegetation_Forestry", "ForestrySites_OP.shp"))
map_WRS_sites <- function(park, title) {
cols <- c("1" = "orangered3", "0" = "wheat", "forest" = "darkseagreen", "grassland" = "lightgoldenrod", "wetland" = "lightblue", "water" = "skyblue4", "LBR" = "red3", "Mature forest" = "darkgreen", "Natural succession" = "blue3", "Seedling" = "purple3")
ggplot() +
geom_sf(data = spatial.boundaries.parks, fill = "grey97", size = 0) +
geom_sf(data = basemap_habs, aes(fill = type_simplified), color = NA, show.legend = FALSE, alpha = .7) +
geom_sf(data = forestry_sites, fill = NA, size = .5) +
geom_sf(data = sites_for_map, aes(color = Type)) +
geom_sf(data = spatial.boundaries.parks, fill = NA, size = .7) +
scale_fill_manual(values = cols) +
scale_fill_manual(values = cols, aesthetics = "color") +
scale_y_continuous(breaks = seq(40, 50, by = .02)) +
scale_x_continuous(breaks = seq(-95, -90, by = .02)) +
theme_minimal() +
theme(
axis.title = element_blank(),
plot.title = element_text(size = 11, face = "bold"),
plot.subtitle = element_text(size = 10, face = "plain")
) +
geom_text_repel(data = sites_for_map %>% filter(Park == ifelse(park == "Cleary", "Cleary Lake", park)), mapping = aes(color = Type, label = sub(".*\\-", "", ExistingID), geometry = geometry), stat = "sf_coordinates", min.segment.length = 1, show.legend = FALSE) +
map_extent_park(park) +
annotation_scale(location = "br", style = "bar", bar_cols = c("grey40", "grey40"), line_col = "grey40", text_col = "grey40", line_width = .5) +
labs( # title = "Sites",
subtitle = paste0(title, "\n", park)
)
}
map_WRS_sites(park = "Lake Rebecca", title = "Woodland Restoration Survey Sites")
map_WRS_sites(park = "Baker", title = "Woodland Restoration Survey Sites")
map_WRS_sites(park = "Crow-Hassan", title = "Woodland Restoration Survey Sites")
map_WRS_sites(park = "Cleary", title = "Woodland Restoration Survey Sites")# Generate distance matrix with default settings (bray-curtis)
dist_mat <- vegdist(wrs_species)
# Determine best dimensionailty by comparing stress levels with scree plot
screeplot <- dimcheckMDS(wrs_species,
distance = "bray", # distance measure used for metaMDS analysis
trymax = 200,
k = 6
) # maximum number of dimensions, default is 6. could in theory test up to max number of species.Stress with just 2 dimensions is >0.2, which is unusable according to some. Going to use 3 dimensions.
set.seed(1)
nmds3 <- metaMDS(
comm = wrs_species, # Define the data to use
distance = "bray", # Specify distance
k = 3,
try = 1000, trymax = 1000, # Number of random starts
autotransform = TRUE, # allow default for community data
noshare = TRUE, # allow default for community data
wascores = TRUE
) # allow default for community data
nmds3## [1] 3
## [1] 11
# nmds3$points # scores for each location (as samples)
# nmds3$species # scores for each animal order (as variables)A Shepard diagram is a scatterplot of distances between points in an nMDS final configuration (observed dissimilarity) against the original dissimilarities (ordination distance). For a nMDS, the line of best fit produced is monotonic, in that it can only step up or remain constant. The points in a good Shepard diagram should be a clean curve or straight line. Large scatter around the line suggests that original dissimilarities are not well preserved in the reduced number of dimensions. If your plot resembles a step-wise or inverse L-shaped function, it would be wise to not interpret this data and instead try to solve for a degenerate solution (not covered in this course). However, please note if you have a small data set (not many observations) your data may look step-wise in a Shepard diagram. Note, you will need to produce a new Shepard diagram when altering the number of dimensions being plotted. Confusingly, the Shepard diagram is called stressplot in R!
TODO: Looks like I might have a step. How problematic is this and what to do about it?
What happens if I run on version without species only observed once?
sp_singletons <- wrs %>%
group_by(Sp) %>%
summarize(n = n()) %>%
filter(n == 1) %>%
pull(Sp)
wrs_species_no_singletons <- wrs %>%
filter(Sp %!in% sp_singletons) %>%
group_by(Survey_ID_full, Sp) %>%
summarize(Count = n()) %>%
spread(Sp, Count, fill = 0) %>%
arrange(Survey_ID_full) %>%
ungroup() %>%
dplyr::select(-Survey_ID_full)
screeplot_no_singletons <- dimcheckMDS(wrs_species_no_singletons,
distance = "bray", # distance measure used for metaMDS analysis
trymax = 200,
k = 6
) # maximum number of dimensions, default is 6. could in theory test up to max number of species.# screeplot_no_singletons
# Stress from just 2 dimensions is still >0.2. going to go with 3 dimensions.
nmds_no_singletons <- metaMDS(
comm = wrs_species_no_singletons, # Define the data to use
distance = "bray", # Specify distance
k = 3,
try = 1000, trymax = 1000, # Number of random starts
autotransform = TRUE, # allow default for community data
noshare = TRUE, # allow default for community data
wascores = TRUE
) # allow default for community data
# nmds_no_singletons
stressplot(nmds_no_singletons)Eliminating singletons didn’t seem to change the Shepard diagram much, so I’m going to use full dataset for rest of analysis.
Basic biplot with ellipses around all sites of each type.
# Define colors for treatments
colvec <- c("green4", "orange", "maroon", "royalblue")
scl <- 3 # scaling term
# Generate biplot with ellipses
plot(nmds3, type = "n", las = 1)
text(nmds3, display = "species", cex = 0.55, col = "grey40")
points(nmds3,
display = "sites", col = colvec,
scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
ordiellipse(nmds3, wrs_surveys$Type,
display = "sites",
kind = "ehull", label = F,
lwd = 2, col = colvec
)
legend("topright",
legend = levels(wrs_surveys$Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec, cex = .5
)Biplot where ellipses represent error around centroid: “When kind=‘sd’ or ‘se’, the ellipse represents a calculated region of error around each group’s centroid. The confidence level is set with conf= and can be understood as how single-dimensional error bars or confidence intervals might be used in univariate plots to estimate the precision of estimated mean or centroid (se) or connect variation in the sampled data to the spread of the broader population (sd)” (Source)
plot(nmds3, type = "n", las = 1)
ordiellipse(nmds3, wrs_surveys$Type,
display = "sites",
kind = "se", conf = 0.95,
label = F, col = colvec,
draw = "polygon", alpha = 50, scaling = scl
)
points(nmds3,
display = "sites", col = colvec,
scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
text(nmds3, display = "species", scaling = scl, cex = 0.55, col = "grey40")
legend("topright",
legend = levels(wrs_surveys$Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec, cex = .5
)Now generate same plot but with filtered species for clarity:
plot(nmds_no_singletons, type = "n", las = 1)
ordiellipse(nmds_no_singletons, wrs_surveys$Type,
display = "sites",
kind = "se", conf = 0.95,
label = F, col = colvec,
draw = "polygon", alpha = 50, scaling = scl
)
points(nmds_no_singletons,
display = "sites", col = colvec,
scaling = scl, pch = 21, cex = 0.75, bg = colvec
)
text(nmds_no_singletons, display = "species", scaling = scl, cex = 0.55, col = "grey40")
legend("topright",
legend = levels(wrs_surveys$Type), bty = "n",
col = colvec, pch = 21, pt.bg = colvec, cex = .5
)Mature forest separates from restored sites (aligned with species like WOTH, EAWP, HAWO, SCTA, OVEN). Restoration types have overlap here, so not major differences (but doesn’t account for site age, etc.).
“The key function for testing groups (and gradients) against an ordination ad hoc–after the ordination has been fit and the latent variable described–is envfit. It is set up similarly to the ordi* functions, but since it is a test, it uses formula notation.”
envfit(nmds3 ~ wrs_surveys$Type + wrs_surveys$Time, choices = c(1:3)) # use all 3 dimensions when hypothesis testing (default is just 2)##
## ***VECTORS
##
## NMDS1 NMDS2 NMDS3 r2 Pr(>r)
## wrs_surveys$Time -0.851510 -0.038286 0.522940 0.4408 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
##
## ***FACTORS:
##
## Centroids:
## NMDS1 NMDS2 NMDS3
## wrs_surveys$TypeMature forest -0.4078 -0.0574 0.1531
## wrs_surveys$TypeLBR 0.0404 -0.1170 -0.2384
## wrs_surveys$TypeSeedling 0.2338 -0.0248 -0.0442
## wrs_surveys$TypeNatural succession 0.1155 0.1502 0.0077
##
## Goodness of fit:
## r2 Pr(>r)
## wrs_surveys$Type 0.2257 0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999
Test indicates that the avian community is significantly different in at least 1 of the woodland types. Time since revegetation also important. Note that envfit tests each factor independently.
We now check pairwise comparisons to see which woodland types differ.
##
## Pairwise comparisons using factor fitting to an ordination
##
## data: nmds3 by wrs_surveys$Type
## 999 permutations
##
## Mature forest LBR Seedling
## LBR 0.004 - -
## Seedling 0.003 0.170 -
## Natural succession 0.003 0.092 0.092
##
## P value adjustment method: fdr
Mature forest differs significantly from all other types. LBR and seedling not significantly different. Some suggestion that both LBR and Seedling sites differ from Natural Succession sites, but not significant at p < 0.05. (p = 0.084).
It would be good to analyze interaction between time and type. I think adonis can be used for this.
Type is highly significant (communities in at least one forest type are different from others). Time alone isn’t significant after controlling for type, but there is an interaction effect (influence of time is mediated by type). Suggests different trajectories for different types…
I don’t fully understand this command, but by = “onedf” has interesting output. From help files “by =”terms” will assess significance for each term (sequentially from first to last), setting by = “margin” will assess the marginal effects of the terms (each marginal term analysed in a model with all other variables), by = “onedf” will analyse one-degree-of-freedom contrasts sequentially, by = NULL will assess the overall significance of all terms together. The arguments is passed on to anova.cca.”
This is kinda interesting. Suggests significant interaction with time and seedling but not LBR (because LBR jump starts process and less variability over time?). Also LBR not signifcantly diff.
Pairwise comparisons:
## $parent_call
## [1] "wrs_species ~ Type , strata = Null , permutations 999"
##
## $`Natural succession_vs_Mature forest`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 1.7598 0.11425 6.9651 0.001 ***
## Residual 54 13.6440 0.88575
## Total 55 15.4039 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Natural succession_vs_Seedling`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.5021 0.03035 1.8778 0.028 *
## Residual 60 16.0431 0.96965
## Total 61 16.5452 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Natural succession_vs_LBR`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.431 0.0412 1.6329 0.041 *
## Residual 38 10.029 0.9588
## Total 39 10.460 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Mature forest_vs_Seedling`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 2.4625 0.13554 9.7208 0.001 ***
## Residual 62 15.7057 0.86446
## Total 63 18.1681 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Mature forest_vs_LBR`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 1.1478 0.10589 4.7372 0.001 ***
## Residual 40 9.6919 0.89411
## Total 41 10.8397 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Seedling_vs_LBR
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.3466 0.02787 1.3187 0.175
## Residual 46 12.0909 0.97213
## Total 47 12.4376 1.00000
##
## attr(,"class")
## [1] "pwadstrata" "list"
## $parent_call
## [1] "wrs_species ~ Type * Time , strata = Null , permutations 999"
##
## $`Natural succession_vs_Mature forest`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 1.7598 0.11425 6.9163 0.001 ***
## Time 1 0.1583 0.01028 0.6221 0.844
## Residual 53 13.4857 0.87548
## Total 55 15.4039 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Natural succession_vs_Seedling`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.5021 0.03035 1.9138 0.027 *
## Time 1 0.2610 0.01578 0.9949 0.473
## Type:Time 1 0.5655 0.03418 2.1555 0.011 *
## Residual 58 15.2166 0.91970
## Total 61 16.5452 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Natural succession_vs_LBR`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.4310 0.04120 1.6146 0.056 .
## Time 1 0.1529 0.01462 0.5730 0.902
## Type:Time 1 0.2678 0.02560 1.0032 0.427
## Residual 36 9.6085 0.91858
## Total 39 10.4602 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Mature forest_vs_Seedling`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 2.4625 0.13554 9.9890 0.001 ***
## Time 1 0.6682 0.03678 2.7107 0.003 **
## Residual 61 15.0375 0.82768
## Total 63 18.1681 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $`Mature forest_vs_LBR`
## Df SumOfSqs R2 F Pr(>F)
## Type 1 1.1478 0.10589 4.7473 0.001 ***
## Time 1 0.2624 0.02421 1.0853 0.367
## Residual 39 9.4295 0.86990
## Total 41 10.8397 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $Seedling_vs_LBR
## Df SumOfSqs R2 F Pr(>F)
## Type 1 0.3466 0.02787 1.3666 0.148
## Time 1 0.5011 0.04029 1.9757 0.021 *
## Type:Time 1 0.4295 0.03453 1.6933 0.045 *
## Residual 44 11.1603 0.89731
## Total 47 12.4376 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## attr(,"class")
## [1] "pwadstrata" "list"
Restoration sites vs reference: Mature forest differs from all other types. Bird communities in restoration sites have not achieved a composition similar to the mature reference sites.
Active vs. passive restoration: Natural succession and seedling sites significantly different. Difference between LBR and natural succession sites is marginally significant. There is a significant interaction between type and time when comparing natural succession and seedling sites (community changes at different rates in seedling sites than natural succession sites), but not for LBR.
Large bareroot vs seedling plantings: No significant difference based on restoration type alone, but the significant interaction between type and time suggests restoration type does effect the rate at which species composition changes.
From adnois help file: “Anderson (2001, Fig. 4) warns that the method may confound location and dispersion effects: significant differences may be caused by different within-group variation (dispersion) instead of different mean values of the groups (see Warton et al. 2012 for a general analysis). However, it seems that adonis2 is less sensitive to dispersion effects than some of its alternatives (anosim, mrpp). Function betadisper is a sister function to adonis2 to study the differences in dispersion within the same geometric framework.”
Also see end of this help vid: https://www.youtube.com/watch?v=oLf0EpMJ4yA
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.02460 0.0081993 0.9149 999 0.416
## Residuals 100 0.89621 0.0089621
Doesn’t seem to be different within-group variation, which is good, I think.
# Tidy distance matrix to use in ggplot
dist_mat_tidy <- dist_mat %>%
as.matrix() %>%
as_tibble(rownames = "samples", ) %>%
pivot_longer(-samples) %>%
left_join(wrs_surveys %>% rownames_to_column("samples"), by = join_by(samples)) %>%
left_join(wrs_surveys %>% rownames_to_column("name"), by = c("name")) %>%
mutate(
comparison = paste(Type.x, Type.y, sep = "-"),
time_difference = Time.x - Time.y
) %>%
rowwise() %>%
mutate(comparison_alpha = paste(sort(c(Type.x, Type.y)), collapse = "-")) %>%
ungroup() %>%
mutate(
samples = as.numeric(samples),
name = as.numeric(name),
combo = paste0(pmin(samples, name), "-", pmax(samples, name))
) %>%
separate(combo, into = c("site_min", "site_max"), sep = "-", remove = FALSE, convert = TRUE) %>%
filter(samples != name) %>% # drop same-site comparisons
filter(time_difference >= 0) %>% # for reciprocal distance measures only keep the positive time distance ()
distinct(site_min, site_max, .keep_all = TRUE) # only keep one when reciprocal sites are different but without time difference, which fel lthrough cracks of previous filters
# filter(samples < name)Each point in these plots is a paired comparison between each site and each mature forest sites. Site pairs with lower Bray-Curtis differences have more similar bird communities (less difference). A values of 1 indicates that the sites have no species in common; a value of 0 indcates that the sites have the same species composition.
# Jitter plot. All woodland types vs mature forest pots.
dist_mat_tidy %>%
filter(str_detect(comparison, "Mature")) %>%
ggplot(aes(x = comparison, y = value)) +
geom_jitter(aes(color = as.numeric(time_difference)), size = .5, width = 0.4, alpha = .9) +
stat_summary(
fun.data = median_hilow, color = "red", size = 1,
fun.args = list(conf.int = 0.50)
) +
labs(x = NULL, y = "Bray-Curtis differences", color = "Time difference between sites (years)") +
# scale_x_discrete(breaks=c("early", "late", "same"),
# labels=c("Inter-mouse\ndistances at\n9 dpw",
# "Inter-mouse\ndistances at\n141 dpw",
# "Intra-mouse\ndistances between\n9 and 141 dpw")) +
scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
theme_classic() +
scale_color_viridis_c() +
theme(axis.text.x = element_text(angle = 90))The mature forest vs. mature forest comparison is a measurement of normal community variability within reference sites (even mature forest sites differ form one another to a certain degree). On average the LBR sites are more similar to mature forests than the seedling and natural succession sites, but the effect is small and significance not tested.
Similar plot but instead of coding time differences with color we plot time difference on the x axis. Time difference is the difference in site ages in years, ie how much older the mature forests are than the restoration).
# Time difference on x-axis. Tall facets.
dist_mat_tidy %>%
filter(str_detect(comparison, "Mature")) %>%
ggplot(aes(x = time_difference, y = value)) +
geom_point(aes(color = Type.y), size = .9, alpha = .25) +
geom_smooth(method = lm, se = TRUE, color = "black") +
stat_summary(
fun.data = median_hilow, color = "grey30", size = .3,
fun.args = list(conf.int = 0.50)
) +
# stat_summary(fun.data=median_hilow, color="red", size=1, fun.args = list(conf.int=0.50)) +
labs(x = NULL, y = "Bray-Curtis differences") +
scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
scale_x_reverse() +
facet_wrap(~Type.y, nrow = 1) +
ggthemes::theme_clean() +
labs(
x = "Time difference (mature forest head-start)",
y = "Community dissimilarity (Bray-Curtis distance)",
color = "Woodland plot type\ncompared to\nmature forest plot"
)
# Time difference on x-axis. Long facets.
dist_mat_tidy %>%
filter(str_detect(comparison, "Mature")) %>%
ggplot(aes(x = time_difference, y = value)) +
geom_point(aes(color = Type.y), size = .9, alpha = .25) +
geom_smooth(method = lm, se = TRUE, color = "black") +
stat_summary(
fun.data = median_hilow, color = "grey30", size = .3,
fun.args = list(conf.int = 0.50)
) +
scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
scale_x_reverse() +
facet_wrap(~Type.y, nrow = 4) +
ggthemes::theme_clean() +
labs(
x = "Time difference (mature forest head-start)",
y = "Plot community dissimilarity (vs. mature forest plots)",
color = "Woodland plot type\ncompared to\nmature forest plot"
)
# Time on x-axis. Overlayed. With "cohorts".
dist_mat_tidy %>%
filter(str_detect(comparison, "Mature")) %>%
ggplot(aes(x = time_difference, y = value)) +
geom_point(aes(color = Type.y), size = .9, alpha = .1) +
geom_smooth(aes(group = Type.y, color = Type.y), method = lm, se = FALSE, size = 2) +
stat_summary(
fun.data = median_hilow, aes(fill = Type.y), color = "grey20", size = .4, pch = 24, stroke = .5,
fun.args = list(conf.int = 0.50)
) +
labs(x = NULL, y = "Bray-Curtis differences") +
scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
scale_x_reverse() +
theme_classic() +
labs(
x = "Time difference (mature forest head-start)",
y = "Plot community dissimilarity (vs. mature forest plots)"
)
# Time on x-axis. Overlayed. Without "cohorts".
dist_mat_tidy %>%
filter(str_detect(comparison, "Mature")) %>%
ggplot(aes(x = time_difference, y = value)) +
geom_point(aes(color = Type.y), size = 1.75, alpha = .1) +
geom_smooth(aes(group = Type.y, color = Type.y), method = lm, se = TRUE, size = 1.25) +
geom_smooth(method = lm, se = FALSE, color = "grey40", size = .5) +
labs(x = NULL, y = "Bray-Curtis differences") +
scale_y_continuous(limits = c(.2, 1), breaks = seq(.2, 1, 0.2)) +
scale_x_reverse() +
theme_classic() +
labs(
x = "Time difference (mature forest head-start)",
y = "Plot community dissimilarity (vs. mature forest plots)",
color = "Woodland plot type\ncompared to\nmature forest plot"
)These plots show that restoration sites and natural succession sites are become more similar to mature forest as they age. LBR sites start out more similar to mature forests than seedling sites. It takes about 20 years for the seedling sites to catch up with LBR sites. Natural succession community progression towards mature sites is slower than LBR and seedling.
By Sam Safran